####This program was used in a BATCH mode
###R CMD BATCH "--args num1 num2 R pp" Table4a.R
## Out-of-sample evaluation, rolling scheme
## It computes the simulations num1 to num2 (for parallelization)
## R is the size of the sample used for estimation
## pp is the risk exposure of the VaR forecast
variable <- commandArgs(trailingOnly=TRUE)
ndeb <- as.numeric(variable[1])
nfin <- as.numeric(variable[2])
r <- as.numeric(variable[3])
pp <- as.numeric(variable[4])


omega0=0.00001
alpha0=0.045
beta0=0.95
student0=6
skew0=0.5
maxp=500

library(fGarch)
library(MASS)
library(nloptr)



#### SImulation Skew T-GARCH
simu_garch_student<-function(nn){
  spec = garchSpec(model = list(omega=omega0, alpha =alpha0, beta = beta0,shape=student0,skew=skew0), cond.dist = "sstd")
  return(garchSim(spec, n = nn))
}

### log p.d.f.
vrstu<-function(x,nu){ 
	lknu=lgamma((nu+1)/2)-0.5*log(pi*(nu-2))-lgamma(nu/2);
	likeli=lknu-((nu+1)/2)*log(1+(x^2)/(nu-2))
return(likeli)
}


### Score Standardized student

scorestu<-function(x,nu){
	knu=0.5*(-1/(nu-2)+digamma((nu+1)/2)-digamma(nu/2));
	fst=(nu+1)/(2*(nu-2))*((x^2)/(nu-2+x^2))-log(1+(x^2)/(nu-2))/2
	return(fst+knu)
}

    ##Log-likelihood
    	garchLLH = function(x,parm) {
        omega = parm[1]; alpha = parm[2]; beta = parm[3];nu=parm[4];
        z =x; ltv = mean(z^2)
        # Use Filter Representation:
        e = omega + alpha * c(ltv, z[-length(x)]^2)
        h = filter(e, beta, "r", init =ltv )
        llh = sum(0.5*log(h)-vrstu(z/sqrt(h),nu))/length(x);
        return(llh)
	}
 

test<-function(ns,r,pp){
  set.seed(10*ns)
  ############## DGP here ###################
    	## We simulate a T- GARCH distribution of size n
	n=r+maxp
    	data=simu_garch_student(n)
    	data=as.vector(data)

	    
	


##### Rolling estimator ######
rolling<-function(h){

    data_est=data[h:(r+h-1)]

    #Log likekihod with constraint  on the parameters

   	
	#Initial conditions
	voldeb=var(data_est)
    	adeb=0.05
    	bdeb=0.9
    	odeb=voldeb*(1-adeb-bdeb);
	
	vraisnew=function(x,b,vol){
		bb=b;
		bb[1]=exp(-b[1]^2)*vol
		bb[3]=0.999*exp(-b[3]^2)
		bb[2]=exp(-b[2]^2)*(0.9999-bb[3])+0.00001
		bb[4]=4.01+96*exp(-b[4]^2)
		return(garchLLH(x,bb))
	}
	

	k4=mean(data_est^4)/(voldeb^2)
	if (k4>3){nudeb=4.011+6/(k4-3)}else
	{nudeb=4.1}
	if (nudeb>100){nudeb=100}
	tdeb=c(sqrt(-log(1-adeb-bdeb)),sqrt(-log((adeb-0.00001)/(0.9999-bdeb))),sqrt(log(0.999/bdeb)),sqrt(-log((nudeb-4.01)/96)))

   	######MLE 
    	log_lik2=function(param){vraisnew(data_est,param,voldeb)}
    	options=list("algorithm"="NLOPT_LN_NELDERMEAD","xtol_rel"=1e-4,"maxeval"=1000)
	res<-nloptr(x0=tdeb,eval_f=log_lik2,opts=options)
	estimpar=res$solution
	omega=voldeb*exp(-estimpar[1]^2)
	beta=0.999*exp(-estimpar[3]^2)
	alpha=exp(-estimpar[2]^2)*(0.9999-beta)+0.00001
    	nu=4.01+96*exp(-estimpar[4]^2)
        nu=round(nu)
        if (nu<=4){nu=4.1}
	
  sauve=c(res$status,omega,alpha,beta,nu)
    ###### Calculation of the score function######
	####dlog sigma/dtheta#########
	#### See Zakoian er al. ######
    	zz=data[h:(r+h)]
	ltv = mean(zz^2)
    	e = omega + alpha * c(ltv, zz[-length(zz)]^2)
    	hh = filter(e, beta, "r", init =ltv )
    	c3=filter(c(0,hh[-(r+1)]),beta,"r",init=0)
    	c2=filter(c(0,zz[-(r+1)]^2),beta,"r",init=0)
    	hh=as.vector(hh)
    	### The derivative of log vol with respect to the parameters.
    	der=cbind(1/hh/(1-beta),c2/hh,c3/hh)
    	der[1,]=c(0,0,0)	

    	######Hit sequence for alpha=p
    	epst=as.vector(as.numeric((zz[1:r]/sqrt(hh[1:r]))))
	   #### Hits
    	It=1*(pt(epst*sqrt(nu/(nu-2)),nu)<pp)
    	VAR=quantile(zz[1:r],pp)
	#It=1*(zz[1:r]<VAR)
	   qa=qt(pp,nu)/sqrt(nu/(nu-2))
	   coef=qt(pp,nu)*dt(qt(pp,nu),nu)

     ### Correction wrt the true score function
	sc=((nu+1)*(epst^2)/(nu-2+epst^2)-1)
    	matV=(2*nu)/(nu+3)
   
        matP=-coef
    	dlog=der[1:r,]/2
 	espder=apply(dlog,2,mean)
 	Ahat=-coef*t(espder)
        score=dlog*sc
    	Vhat=ginv(t(score)%*%score/r)
    	corscore=(Ahat)%*%Vhat%*%t(Ahat)

	##### Oblique projection, tilde(e)_t in the paper
    	etildet=It-pp-matP/2*(epst^2-1)
	   vett=as.numeric(mean(etildet^2))

	############### Corrected hits ###########################
	############## sans nu ##########
    	Itloc=It-pp	
    	Itm1=c(NA,Itloc[-length(It)])
    	Itm2=c(NA,Itm1[-length(It)])
    	Itm3=c(NA,Itm2[-length(It)])
    	Ahat1=c(mean(Itloc*Itm1*score[,1],na.rm=TRUE),mean(Itloc*Itm1*score[,2],na.rm=TRUE),
		  mean(Itloc*Itm1*score[,3],na.rm=TRUE))
    	corscore1=t(Ahat1)%*%Vhat%*%Ahat1

    	Ahat2=c(mean(Itloc*Itm2*score[,1],na.rm=TRUE),mean(Itloc*Itm2*score[,2],na.rm=TRUE),
		  mean(Itloc*Itm2*score[,3],na.rm=TRUE))
    	corscore2=t(Ahat2)%*%Vhat%*%Ahat2

    	Ahat3=c(mean(Itloc*Itm3*score[,1],na.rm=TRUE),mean(Itloc*Itm3*score[,2],na.rm=TRUE),
		  mean(Itloc*Itm3*score[,3],na.rm=TRUE))
    	corscore3=t(Ahat3)%*%Vhat%*%Ahat3

    #VaR forecasts from R+1
    volf=hh[r+1] # Variance forecast
    epstf=zz[r+1]/sqrt(volf) # ratio of real return to vol forecast
    Itf=1*(pt(epstf*sqrt(nu/(nu-2)),nu)<pp)
    #Itf=1*(zz[r+1]<VAR)
	############## Our moments ##########
	scf=((nu+1)*(epstf^2)/(nu-2+epstf^2)-1)
	dlogf=der[(r+1),]/2
    scoref=dlogf*scf

	###### e_t in the paper, i.e. projection onto the orthogonal of an auxiliary score
    	et=Itf-pp+coef*scf/matV
    	# Normalization to have unit variance
    	etf=et/as.numeric(sqrt(pp*(1-pp)-coef^2/matV))

    	##### etstar in the paper, i.e. projection onto the orthogonal of the true score
    	etstar=Itf-pp-scoref%*%Vhat%*%t(Ahat)
    	# Normalization to have unit variance
    	etstf=etstar/as.numeric(sqrt((pp*(1-pp)-corscore)))

    	##### Correction of the variance of the hit sequence
    	Itbis=(Itf-pp)/as.numeric(sqrt((pp*(1-pp)-corscore)))

    	##### Oblique projection, tilde(e)_t in the paper
    	etildet=Itf-pp-matP[1]/2*(epstf^2-1)
	   etfbis=etildet/sqrt(vett)

	weight<-function(p,r){
		if ((p/r)<1){poids1=-(p/r)^2/3}else{poids1=(2*r)/3/p-1}
		return(poids1)
	}
	calculcov<-function(pp,cor,p,r){
		poids=weight(p,r)
    		return((pp*(1-pp)*pp*(1-pp)+poids*cor)^{0.25})
	}
	calculv<-function(pp,cor,p,r){
		poids=weight(p,r)
    		return((pp*(1-pp)+poids*cor)^{0.5})
	}

	Itcf=Itf-pp
	return(cbind(etf,etfbis,etstf,Itcf,Itcf/calculv(pp,corscore,100,r),Itcf/calculv(pp,corscore,250,r)
,Itcf/calculv(pp,corscore,500,r),Itcf/calculcov(pp,corscore1,100,r),Itcf/calculcov(pp,corscore1,250,r)
,Itcf/calculcov(pp,corscore1,500,r),Itcf/calculcov(pp,corscore2,100,r),Itcf/calculcov(pp,corscore2,250,r)
,Itcf/calculcov(pp,corscore2,500,r),Itcf/calculcov(pp,corscore3,100,r),Itcf/calculcov(pp,corscore3,250,r)
,Itcf/calculcov(pp,corscore3,500,r)))
}
	matrol=t(sapply(seq(1,maxp,1),rolling))

	
	calcultest<-function(vec,nbp){
		x=vec[1:nbp]
		xm1=c(NA,x[-length(x)])
    		xm2=c(NA,xm1[-length(x)])
    		xm3=c(NA,xm2[-length(x)])
      	xix1=(nbp-1)*mean(xm1*x,na.rm=TRUE)^2
    		xix2=(nbp-2)*mean(xm2*x,na.rm=TRUE)^2
    		xix3=(nbp-3)*mean(xm3*x,na.rm=TRUE)^2
		xiM=nbp*mean(x)^2
		return(c(xiM,xix1,xix2,xix3))
	}
	##### test P=100
	testet=calcultest(matrol[,1],100)
	testetilde=calcultest(matrol[,2],100)
    	testetstar=calcultest(matrol[,3],100)
	testIt=100*mean(matrol[1:100,5])^2
	
	x=matrol[1:100,8]
	xm1=c(NA,x[-length(x)])
	testItm1=99*mean(xm1*x,na.rm=TRUE)^2

	x=matrol[1:100,11]
	xm1=c(NA,x[-length(x)])
	xm2=c(NA,xm1[-length(x)])
	testItm2=98*mean(xm2*x,na.rm=TRUE)^2

	x=matrol[1:100,14]
	xm1=c(NA,x[-length(x)])
	xm2=c(NA,xm1[-length(x)])
	xm3=c(NA,xm2[-length(x)])
	testItm3=97*mean(xm3*x,na.rm=TRUE)^2

	test100=c(testet[1],testetilde[1],testetstar[1],testIt,testet[2:4],testetilde[2:4],testetstar[2:4],testItm1,testItm2,testItm3)
	##
	##### test P=250
	testet=calcultest(matrol[,1],250)
	testetilde=calcultest(matrol[,2],250)
    	testetstar=calcultest(matrol[,3],250)
	testIt=250*mean(matrol[1:250,6])^2
	
	x=matrol[1:250,9]
	xm1=c(NA,x[-length(x)])
	testItm1=249*mean(xm1*x,na.rm=TRUE)^2

	x=matrol[1:250,12]
	xm1=c(NA,x[-length(x)])
	xm2=c(NA,xm1[-length(x)])
	testItm2=248*mean(xm2*x,na.rm=TRUE)^2

	x=matrol[1:250,15]
	xm1=c(NA,x[-length(x)])
	xm2=c(NA,xm1[-length(x)])
	xm3=c(NA,xm2[-length(x)])
	testItm3=247*mean(xm3*x,na.rm=TRUE)^2

	test250=c(testet[1],testetilde[1],testetstar[1],testIt,testet[2:4],testetilde[2:4],testetstar[2:4],testItm1,testItm2,testItm3)
	

	##### test P=500
	testet=calcultest(matrol[,1],500)
	testetilde=calcultest(matrol[,2],500)
    	testetstar=calcultest(matrol[,3],500)
	testIt=500*mean(matrol[1:500,7])^2
	
	x=matrol[1:500,10]
	xm1=c(NA,x[-length(x)])
	testItm1=499*mean(xm1*x,na.rm=TRUE)^2

	x=matrol[1:500,13]
	xm1=c(NA,x[-length(x)])
	xm2=c(NA,xm1[-length(x)])
	testItm2=498*mean(xm2*x,na.rm=TRUE)^2

	x=matrol[1:500,16]
	xm1=c(NA,x[-length(x)])
	xm2=c(NA,xm1[-length(x)])
	xm3=c(NA,xm2[-length(x)])
	testItm3=497*mean(xm3*x,na.rm=TRUE)^2

	test500=c(testet[1],testetilde[1],testetstar[1],testIt,testet[2:4],testetilde[2:4],testetstar[2:4],testItm1,testItm2,testItm3)
	    return(c(ns,r,pp,100,test100,250,test250,500,test500))
    }


### We launch the program and display the results
options(width=1000)
toto=(matrix(sapply(ndeb:nfin,test,r=r,p=pp),nrow=54))
noquote(formatC(t(toto),digits=3,format="f"))

